---
title: "The Broader Political Significance of Houses of Worship: Online Appendix"
author: |
  | Nathanael Gratias Sumaktoyo
  | National University of Singapore
  | nathanael.sumaktoyo@fulbrightmail.org
output:
  officedown::rdocx_document:
  reference_docx: online_appendix.dotx
toc: yes
toc_depth: 2
number_sections: yes
df_print: paged
---
  
  
version: `r format(Sys.time(), '%d %B, %Y -- %H:%M:%S')`

```{r, echo=FALSE, warning=FALSE, results='hide', include=FALSE}
rm(list=ls())
source("extras/packages.R")
source("extras/functions.R")
library("flextable")
# webshot::install_phantomjs(force = TRUE)

knitr::opts_chunk$set(rows.print=100)

```


\newpage

# Number of New Mosques, by District

<br>

The number of new mosques varies significantly across the treated districts, ranging from 1 to 54. The analysis takes this variation into account through a robustness check that defines treatment as having exactly one new mosque. Please refer to the main text.

<br>
  
```{r, echo=FALSE, warning = FALSE}

load("data/dat_kec_long.rData")

dat = subset(dat, 
             subset = (religion_07 == 1) & (religion_14 == 1) & 
               (census10_prop_muslim >= .5) & (mig_status == "stayer"))

cols = c("year", "change_mo_07_13", "n_mo_07",
         "nmprov_mosque", "nmkab_mosque", "nmkec_mosque", "kecid")

dat = dat[complete.cases(dat[,cols]), cols]

df = dat %>%
  group_by(kecid) %>%
  dplyr::summarise(change_mo_07_13 = mean(change_mo_07_13, na.rm=TRUE),
                   n_mo_07 = mean(n_mo_07, na.rm=TRUE))

tab = data.frame(table(df$change_mo_07_13))

tab$Var1 = as.character(tab$Var1)
tab$Var1 = as.numeric(tab$Var1)


p = ggplot(tab, aes(x = Var1, y = Freq)) +
  geom_bar(stat = "identity") +
  scale_x_continuous(name = "Number of New Mosques (2007-2013)", 
                     breaks = seq(0, max(tab$Var1), 5)) +
  scale_y_continuous(name = "Number of Districts", breaks = seq(0, 200, 20)) +
  labs(subtitle = "Histogram of New Mosques (2007-2013) in Analyzed Districts") + 
  theme_bw()

p

tab1 = data.frame(table(df$change_mo_07_13))
colnames(tab1) = c("Number of New Mosques", "Frequency")
tab2 = data.frame(round(prop.table(table(df$change_mo_07_13)),3))
colnames(tab2) = c("Number of New Mosques", "Proportion")

tab = cbind(tab1, tab2$Proportion)
colnames(tab) = c("Number of New Mosques", "Frequency", "Proportion")

out <- flextable(tab)
out <- fontsize(out, size = 7, part = "all")
out <- set_table_properties(out, width = .8, layout = "autofit")
out
```


\newpage

# New Mosques by Existing Mosques

<br>

Districts with more mosques prior to 2007 also tended to have a higher number of new mosques built between 2007 and 2013. The analysis takes this into account by including the number of mosques in 2007 as a covariate in robustness checks. Please refer to the main text.

<br>
  
```{r, echo=FALSE}

r = cor(df$change_mo_07_13, df$n_mo_07)


### correlation between existing mosques and new mosques

p = ggplot(df, aes(x = n_mo_07, y = change_mo_07_13)) +
  geom_point(alpha = .4) + 
  scale_x_continuous(name = "Number of Mosques in 2007", 
                     breaks = seq(0, max(df$n_mo_07), 50)) + 
  scale_y_continuous(name = "Number of New Mosques (2007-2013)", 
                     breaks = seq(0, max(df$change_mo_07_13), 5)) +
  labs(subtitle = paste0("Correlation = ", round(r, 3))) + 
  theme_bw()

p

```


\newpage

# Geographic Distribution of Respondents across Kabupatens

<br>

The figure presents the distribution of analyzed respondents across kabupatens (the second-level administrative subdivision just above kecamatans or districts). Visualizing on the kabupaten level is driven by practical considerations as a visualization on the district level is almost illegible due to the very small sizes of the territories. The figure suggests the data offers a good coverage of the country's population centers. Most of the respondents came from the Java and Nusa Tenggara islands, where 55\% of Indonesians lived. In addition, considerable numbers of respondents were also sampled from most of Sumatera, the eastern part of Kalimantan, and the southern part of Sulawesi.

<br>

```{r, echo=FALSE, warning=FALSE, results='hide', message=FALSE, fig.width=6}

rm(list = setdiff(ls(), lsf.str()))

load("data/dat_kec_wide.rData")

dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
               (census10_prop_muslim >= .5) & (mig_status == "stayer"))

cols = c("change_mo_07_13",
         "nmprov_mosque", "nmkab_mosque", "pidlink")

dat = dat[complete.cases(dat[,cols]), cols]
dat = dat[!duplicated(dat$pidlink),]
dat$kabid = paste0(dat$nmprov_mosque, " - ", dat$nmkab_mosque)

dat = dat %>% 
  group_by(nmprov_mosque, nmkab_mosque) %>%
  dplyr::summarise(n = n())

### geospatial data
path = getwd()
map = gadm(country = "IDN", path = path, level = 2)

# to lower case
map$NAME_1 = tolower(map$NAME_1)
map$NAME_2 = tolower(map$NAME_2)

# create nmprov nmkab
map$nmprov_gadm = map$NAME_1
map$nmkab_gadm = map$NAME_2

cw = read.xlsx2("data/crosswalk_kab_gadm_mosque.xlsx", sheetIndex = 1)

### merge gadm and ifls
map = merge(map, cw,
             by = c("nmprov_gadm", "nmkab_gadm"),
                all.x = TRUE, all.y = FALSE)

map = merge(map, dat,
           by = c("nmprov_mosque", "nmkab_mosque"),
                all.x = TRUE, all.y = FALSE)


p = ggplot(map) + 
  geom_spatvector(aes(fill = n)) +
  scale_fill_distiller(name = "", 
                       palette = "YlGn",
                       breaks = pretty_breaks(n = 10),
                       trans = "reverse") +
  labs(subtitle = "Distribution of Analyzed Respondents across Kabupatens")


p

```

\newpage

# Descriptive Statistics of Variables

<br>

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))

library(summarytools)

```

## Pre-Treatment (2007)

<br>


```{r, echo=FALSE}

cols = c("live_village2pt", "live_neighborhood2pt", "live_room2pt", 
         "worship_house2pt", "intermarriage2pt", 
         "similar_religion", "similar_ethnicity", 
         "trust_religion2pt", "trust_ethnicity2pt", "helping2pt", 
         "female", "age", "edu", "how_religious", "shalat3pt", "work_lastyear", "marital", 
         "javaisland", "census10_median_age", "census10_prop_male", "census10_prop_working", 
         "census10_median_edu",
         "census10_prop_muslim", "census10_diversity_lang", 
         "elec14_kec_jokowi_share", "elec19_islamist_share", "elec19_jokowi_share",
         "n_mo_07")

df = dfSummary(dat[dat$year == "2007", cols],
          labels.col = FALSE,
          graph.col = FALSE,
          valid.col = FALSE,
          na.col = FALSE,
          max.distinct.values = 4,
          max.string.width = 30,
          headings = FALSE)

df = data.frame(df)
df = lapply(df, function(x) gsub("\\\\","", x))
df = data.frame(df)

out = flextable(data.frame(df))
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = 1, layout = "autofit")

```

## Post-Treatment (2014)

<br>

```{r, echo=FALSE}

cols = c("live_village2pt", "live_neighborhood2pt", "live_room2pt", 
         "worship_house2pt", "intermarriage2pt", 
         "similar_religion", "similar_ethnicity", 
         "trust_religion2pt", "trust_ethnicity2pt", "helping2pt", 
         "female", "age", "edu", "how_religious", "shalat3pt", "work_lastyear", "marital", 
         "javaisland", "census10_median_age", "census10_prop_male", "census10_prop_working", 
         "census10_median_edu",
         "census10_prop_muslim", "census10_diversity_lang", 
         "elec14_kec_jokowi_share", "elec19_islamist_share", "elec19_jokowi_share",
         "n_mo_14")

df = dfSummary(dat[dat$year == "2014", cols],
          labels.col = FALSE,
          graph.col = FALSE,
          valid.col = FALSE,
          na.col = FALSE,
          max.distinct.values = 4,
          max.string.width = 30,
          headings = FALSE)

df = data.frame(df)
df = lapply(df, function(x) gsub("\\\\","", x))
df = data.frame(df)


out = flextable(data.frame(df))
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = 1, layout = "autofit")
```


\newpage

# Population and Sample Sizes of Analyzed Districts

<br>

More than 50% of the analyzed districts had fewer than 100,000 residents, as per the 2010 Census. About 40% may even be considered small communities as they had fewer than 50,000 residents.

<br>

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")

### kecamatan data
load("data/dat_kec_wide.rData")

dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
               (census10_prop_muslim >= .5) & (mig_status == "stayer"))

cols = c("change_mo_07_13", "kecid", "pidlink")

dat = dat[complete.cases(dat[,cols]), ]

```


```{r, echo=FALSE, warning=FALSE}

df = dat %>% group_by(kecid) %>%
      dplyr::summarise(tot_pop = mean(census10_total_pop, na.rm=TRUE),
                       n_respondent = n(),
                       change_mo_07_13 = mean(change_mo_07_13))

#### total population
tab = data.frame(table(df$tot_pop))
tab$group = NA
tab$group[tab$Var1 %in% 1:50000] = "<= 50,000"
tab$group[tab$Var1 %in% 50001:100000] = "50,001 -\n100,000"
tab$group[tab$Var1 %in% 100001:150000] = "100,001 -\n150,000"
tab$group[tab$Var1 %in% 150001:200000] = "150,001 -\n200,000"
tab$group[tab$Var1 %in% 200001:250000] = "200,001 -\n250,000"
tab$group[tab$Var1 %in% 250001:300000] = "251,001 -\n300,000"
tab$group[tab$Var1 %in% 300001:350000] = "300,001 -\n350,000"
tab$group[tab$Var1 %in% 350001:400000] = "350,001 -\n400,000"
tab$group[tab$Var1 %in% 400001:450000] = "400,001 -\n450,000"
tab$group[tab$Var1 %in% 450001:500000] = "450,001 -\n500,000"
tab$group[tab$Var1 %in% 500001:600000] = "> 500,000"

tab$group = factor(tab$group,
                   levels = unique(tab$group))

tab = tab %>% group_by(group) %>% 
        dplyr::summarise(n = n())

p <- ggplot(tab, aes(x=group, y = n)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = n), hjust = -.2) +
  scale_x_discrete(name = "Population Size (Census 2010)") +
  scale_y_continuous(name = "Number of Analyzed Districts",
                     limits = c(0, 500)) +
  labs(subtitle = "Population Size of Analyzed Districts") +
  theme_bw() +
  coord_flip()

p

tab = summary(df$tot_pop)
tab <- setNames(tab, names(tab))
tab <- as.data.frame(t(tab))[,c(2,3)]

colnames(tab) = c("Statistic", "Value")

out <- flextable(tab)
out <- set_table_properties(out, width = 1, layout = "autofit")
out

```




```{r, echo=FALSE}

#### sample size

tab = data.frame(table(df$n_respondent))
tab$group = NA

start = seq(1, 211, 10)
end = seq(10, 220, 10)

for (i in 1:length(start)) {
  tab$group[tab$Var1 %in% start[i]:end[i]] = paste0(start[i], " - ", end[i])
}

tab$group = factor(tab$group,
                   levels = unique(tab$group))

tab = tab %>% group_by(group) %>% 
        dplyr::summarise(n = sum(Freq))

p <- ggplot(tab, aes(x=group, y = n)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = n), hjust = -.2) +
  scale_x_discrete(name = "Sample Size") +
  scale_y_continuous(name = "Number of Analyzed Districts",
                     limits = c(0,850)) +
  labs(subtitle = "Sample Size in the Analyzed Districts") +
  theme_bw() + coord_flip()

p


tab = summary(df$n_respondent)
tab <- setNames(tab, names(tab))
tab <- as.data.frame(t(tab))[,c(2,3)]

colnames(tab) = c("Statistic", "Value")

out <- flextable(tab)
out <- set_table_properties(out, width = 1, layout = "autofit")
out

```



\newpage

<!---BLOCK_LANDSCAPE_START--->

# Main Regression Models

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))


dat$treat = ifelse(dat$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))


```

<br>

Number of districts by treatment status:

```{r, echo=FALSE}

d = dat[!duplicated(dat[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(d$treat))
tab2 = data.frame(round(prop.table(table(d$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status:

```{r, echo=FALSE}

tab1 = data.frame(table(dat$treat) / 2)
tab2 = data.frame(round(prop.table(table(dat$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")
```


<br>


```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA,
                   n_obs = NA,
                   n_kec = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx],
                   n_obs = nobs(out$m)/2,
                   n_kec = ngrps(out$m)))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  
  
  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_main_models.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")


```

\newpage

# Magnitudes of Effects among Attendees and Non-Attendees

<br>

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))




```

```{r, echo=FALSE}
dat$treat = ifelse(dat$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```


## Interaction: Gender (Female)

<br>

```{r, echo=TRUE}
dat$mod = dat$female
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_ind_cov/female.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```

## Interaction: Proportion of Males in District

<br>

```{r, echo=TRUE}
dat$mod = dat$census10_prop_male
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```

## Interaction: Attending Communal Prayer

<br>

```{r, echo=TRUE}
dat$mod = dat$communal_pray
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
    
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_ind_cov/communalpray.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```

\newpage

# Alternative Specifications

<br>

## Alternative 1 (Binary DV with Covariates)

<br>

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))

dat$treat = ifelse(dat$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```


```{r, echo=FALSE, warning=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  
  
  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_specification_1.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```


## Alternative 2 (Ordinal DV with no Covariates)

<br>


```{r, echo=FALSE, warning=FALSE}

vec_dv = c("live_village", "live_neighborhood", "live_room",
           "worship_house", "intermarriage", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion", "trust_ethnicity", "helping")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_specification_2.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```


## Alternative 3 (Ordinal DV with Covariates)

<br>

```{r, echo=FALSE, warning=FALSE}

vec_dv = c("live_village", "live_neighborhood", "live_room",
           "worship_house", "intermarriage", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion", "trust_ethnicity", "helping")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_specification_3.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```


## Alternative 4 (Districts with n > 10)

<br>

```{r, echo=FALSE}

d = dat[which(!duplicated(dat$pidlink)), ]
d = d %>% group_by(kecid) %>%
          dplyr::summarise(n = n())

d = subset(d, n > 10)

```

Number of districts by treatment status:

```{r, echo=FALSE}

df = dat[dat$kecid %in% d$kecid, ]
df = df[!duplicated(df[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(df$treat))
tab2 = data.frame(round(prop.table(table(df$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status:

```{r, echo=FALSE}

tab1 = data.frame(table(dat[dat$kecid %in% d$kecid,"treat"]) / 2)
tab2 = data.frame(round(prop.table(table(dat[dat$kecid %in% d$kecid,"treat"])),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")
```


<br>

```{r, echo=FALSE, warning=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat[which(dat$kecid %in% d$kecid), ])

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_specification_4.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")


```



## Alternative 5 (Zero or One New Mosque)

<br>

```{r, echo=FALSE}
dat$treat = NA
dat$treat[dat$change_mo_07_13 == 0] = "No New Mosque"
dat$treat[dat$change_mo_07_13 == 1] = "One New Mosque"

dat$treat = factor(dat$treat, levels = c("No New Mosque", "One New Mosque"))
```

Number of districts by treatment status:

```{r, echo=FALSE}

d = dat[!duplicated(dat[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(d$treat))
tab2 = data.frame(round(prop.table(table(d$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status:

```{r, echo=FALSE}

tab1 = data.frame(table(dat$treat) / 2)
tab2 = data.frame(round(prop.table(table(dat$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")


out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")
```


<br>

```{r, echo=FALSE, warning=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatOne New Mosque")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_specification_5.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```

\newpage

<!---BLOCK_LANDSCAPE_STOP--->

## Alternative 6 (Preprocessing with Matching)

<br>

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))


dat$shalat3pt_1 = as.numeric(dat$shalat3pt == 1)
dat$shalat3pt_2 = as.numeric(dat$shalat3pt == 2)
dat$shalat3pt_3 = as.numeric(dat$shalat3pt == 3)


```


```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6}

df = dat[dat$year == "2007", ] %>%
        group_by(kecid) %>%
        dplyr::summarise(
                  javaisland = mean(javaisland, na.rm=TRUE),
                  census10_median_age = mean(census10_median_age, na.rm=TRUE),
                  census10_prop_male = mean(census10_prop_male, na.rm=TRUE),
                  census10_prop_working = mean(census10_prop_working, na.rm=TRUE),
                  census10_median_edu = mean(census10_median_edu, na.rm=TRUE),
                  census10_prop_muslim = mean(census10_prop_muslim, na.rm=TRUE),
                  census10_diversity_lang = mean(census10_diversity_lang, na.rm=TRUE),
                  elec14_kec_jokowi_share = mean(elec14_kec_jokowi_share, na.rm=TRUE),
                  elec19_islamist_share = mean(elec19_islamist_share, na.rm=TRUE),
                  elec19_jokowi_share = mean(elec19_jokowi_share, na.rm=TRUE),
                  ln_n_mo_07 = mean(ln_n_mo_07, na.rm=TRUE),
                  female =  mean(female, na.rm=TRUE),
                  age =  mean(age, na.rm=TRUE),
                  edu =  mean(edu, na.rm=TRUE),
                  how_religious = mean(how_religious, na.rm=TRUE),
                  shalat3pt_1 = mean(shalat3pt_1, na.rm=TRUE),
                  shalat3pt_2 =  mean(shalat3pt_2, na.rm=TRUE),
                  shalat3pt_3 =  mean(shalat3pt_3, na.rm=TRUE),
                  work_lastyear =  mean(work_lastyear, na.rm=TRUE),
                  marital =  mean(marital, na.rm=TRUE),

                  change_mo_07_13 = mean(change_mo_07_13, na.rm=TRUE),

                  live_village2pt = mean(live_village2pt,na.rm=TRUE),
                  live_neighborhood2pt = mean(live_neighborhood2pt,na.rm=TRUE),
                  live_room2pt = mean(live_room2pt,na.rm=TRUE),
                  worship_house2pt = mean(worship_house2pt,na.rm=TRUE),
                  intermarriage2pt = mean(intermarriage2pt,na.rm=TRUE),
                  similar_religion = mean(similar_religion, na.rm=TRUE),
                  similar_ethnicity = mean(similar_ethnicity, na.rm=TRUE),
                  trust_religion2pt = mean(trust_religion2pt,na.rm=TRUE),
                  trust_ethnicity2pt = mean(trust_ethnicity2pt,na.rm=TRUE),
                  helping2pt = mean(helping2pt,na.rm=TRUE))

id_tr0 = which(df$change_mo_07_13 == 0)
id_tr1 = which(df$change_mo_07_13 == 1)

df$treat = NA
df$treat[id_tr0] = 0
df$treat[id_tr1] = 1

df = df[complete.cases(df),]

bal_matched <- matchit(treat ~ 
                       javaisland + census10_median_age + census10_prop_male + 
                         census10_prop_working + census10_median_edu + census10_prop_muslim +
                    census10_diversity_lang + 
                    elec14_kec_jokowi_share + elec19_islamist_share + elec19_jokowi_share +
                      ln_n_mo_07 + 
                    female + age + edu +
                    how_religious + shalat3pt_1 + shalat3pt_2 + shalat3pt_3 + 
                      work_lastyear + marital + 
                    live_village2pt + live_neighborhood2pt + live_room2pt + 
                    worship_house2pt + intermarriage2pt + 
                    similar_religion + similar_ethnicity + 
                    trust_religion2pt + trust_ethnicity2pt + helping2pt,
  data = df, 
  method = "nearest", distance ="glm", link = "probit", discard = "both",
  #ratio = 1,
  replace = FALSE)



out = data.frame(summary(bal_matched)$sum.matched)
out = round(out, 3)
out$Variable = rownames(out)


out$Variable = change_coefs(out$Variable)

out = out[,c("Variable", colnames(out)[1:3])]

out <- flextable(data.frame(out))
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```


<br>


```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6, fig.height=6}

out = data.frame(summary(bal_matched)$nn)
out$Variable = rownames(out)
out = out[,c("Variable", colnames(out)[1:2])]


out <- flextable(data.frame(out))
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

plot(bal_matched, type = "jitter", interactive = FALSE)


idx = which(bal_matched$weights == 1)

df_matched = data.frame(kecid = df$kecid[idx],
                        distance = bal_matched$distance[idx])

```



<!---BLOCK_LANDSCAPE_START--->

```{r, echo=FALSE}
### merge data with propensity score
dat = merge(x = dat, y = df_matched,
            by = "kecid",
            all = FALSE)

dat$treat = ifelse(dat$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))
```

<br>

Number of districts by treatment status:

```{r, echo=FALSE}

d = dat[!duplicated(dat[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(d$treat))
tab2 = data.frame(round(prop.table(table(d$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status: 

```{r, echo=FALSE}
tab1 = data.frame(table(dat$treat) / 2)
tab2 = data.frame(round(prop.table(table(dat$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")


out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

```{r, echo=FALSE, warning=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_matched(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_specification_6.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")


```

\newpage

# Placebo Tests

<br>

## Placebo 1 (In-Space Placebo)


```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")

### kecamatan data
load("data/dat_kec_long.rData")

dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (change_mo_07_13 == 0) &
                            (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))




```

```{r, echo=FALSE}

set.seed(10)
kecids = unique(dat$kecid)
t = sample(kecids, size = length(kecids)/2)

dat$treat = ifelse(dat$kecid %in% t, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))
```

<br>

Number of districts by treatment status:

```{r, echo=FALSE}

d = dat[!duplicated(dat[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(d$treat))
tab2 = data.frame(round(prop.table(table(d$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status: 

```{r, echo=FALSE}
tab1 = data.frame(table(dat$treat) / 2)
tab2 = data.frame(round(prop.table(table(dat$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")


out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

```{r, echo=FALSE, warning=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_placebo_inspace.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```


## Placebo 2 (In-Time Placebo)

<br>


```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")

### kecamatan data
load("data/dat_kec_long.rData")

dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                           (census10_prop_muslim >= .5) & (mig_status == "stayer"))


dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))

dat$treat = ifelse(dat$change_mo_00_06 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```


Number of districts by treatment status:

```{r, echo=FALSE}

d = dat[!duplicated(dat[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(d$treat))
tab2 = data.frame(round(prop.table(table(d$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status:

```{r, echo=FALSE}


tab1 = data.frame(table(dat$treat) / 2)
tab2 = data.frame(round(prop.table(table(dat$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

```{r, echo=FALSE, warning=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_placebo_intime.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```


## Placebo 3 (Analysis of Movers)

<br>


```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                           (census10_prop_muslim >= .5) & (mig_status == "mover"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))



dat$treat = ifelse(dat$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```


Number of districts by treatment status:

```{r, echo=FALSE}

d = dat[!duplicated(dat[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(d$treat))
tab2 = data.frame(round(prop.table(table(d$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status:

```{r, echo=FALSE}

tab1 = data.frame(table(dat$treat) / 2)
tab2 = data.frame(round(prop.table(table(dat$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")


out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

```{r, echo=FALSE, warning=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_placebo_movers.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```


\newpage

# Mechanism: Religiosity 

<br>

## Salat and Religiosity Self-Rating

<br>

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))


dat$treat = ifelse(dat$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```



```{r, echo=FALSE, warning=FALSE}

vec_dv = c("shalat2pt", "how_religious2pt")

vec_dvlab = c("Shalat Five Times or More",
              "Self-Rated Religiosity")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  
 
  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
 
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_religiosity.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "all")
set_table_properties(tab, width = .6, layout = "autofit")

```

## Effects of Mushalla

<br>

```{r, echo=FALSE, warning=FALSE}

rm(list = setdiff(ls(), lsf.str()))

### kecamatan data
load("data/dat_kec_long.rData")

dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))

dat$treat = ifelse(dat$change_mu_07_13 > 0, "New Mushalla(s)", "No New Mushalla")
dat$treat = factor(dat$treat, levels = c("No New Mushalla", "New Mushalla(s)"))

```


Number of districts by treatment status:

```{r, echo=FALSE}

d = dat[!duplicated(dat[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(d$treat))
tab2 = data.frame(round(prop.table(table(d$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status:

```{r, echo=FALSE}

tab1 = data.frame(table(dat$treat) / 2)
tab2 = data.frame(round(prop.table(table(dat$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")


out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")
```

<br>

```{r, echo=FALSE, warning=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mushalla(s)")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  
 
  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
 
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_mushallas.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```


\newpage

# Mechanism: Psychosocial

<br>

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))




```

```{r, echo=FALSE}
dat$treat = ifelse(dat$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```


## Interaction: Friday Interview

<br>

```{r, echo=TRUE}
dat$mod = dat$isfriday
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_ind_cov/friday.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```


## Interaction: NU/Muhammadiyah Affiliation

<br>

```{r, echo=TRUE}
dat$mod = dat$reltrad_numuha
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):modNU/Muhammadiyah")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_ind_cov/numuha.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```



## Interaction: Proportion of Muslims in the District

<br>

```{r, echo=TRUE}
dat$mod = dat$census10_prop_muslim
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_kec_cov/propmuslim.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```



## Interaction: Proportion Waqf Mosques

<br>

```{r, echo=TRUE}
dat$mod = dat$pr_wakaf13_mo
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_kec_cov/propwakaf.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```



## Interaction: Log(Number of Mosques in 2007)

<br>

```{r, echo=TRUE}
dat$mod = dat$ln_n_mo_07
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_kec_cov/nmo07.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```



## Interaction: Proportion NU and Muhammadiyah in District (2014)

<br>

```{r, echo=TRUE}
dat$mod = dat$prop_kec_numuha
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_kec_cov/numuha.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```

\newpage

# Mechanism: Informational

<br>

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))




```

```{r, echo=FALSE}
dat$treat = ifelse(dat$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```


## Interaction: Individual Education

<br>


```{r, echo=TRUE}
dat$mod = dat$edu
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_ind_cov/edu.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```



## Interaction: Cognitive Measure

<br>

```{r, echo=TRUE}
dat$mod = dat$w_abilrc_14

```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_ind_cov/cognitive.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```


## Interaction: Alternative Information Sources


<br>

```{r, echo=TRUE}

dat$papers = as.numeric(dat$papers_idlang | dat$papers_othlang)
dat$mod = rowSums(dat[,c("own_tv_2014", "cellphone_internet_14", "cellphone_socmed_14",
                         "papers")], na.rm=TRUE)

```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_ind_cov/infosources.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```

## Interaction: Median of Education Level (District)

<br>


```{r, echo=TRUE}
dat$mod = dat$census10_median_edu
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_int_kec_cov/medianedu.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```

\newpage

# Other Models Referenced in the Main Text

<br>

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))




```

```{r, echo=FALSE}
dat$treat = ifelse(dat$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```

## Interaction: NU and Muhammadiyah Affiliations

<br>

```{r, echo=TRUE}
dat$mod = dat$reltrad_3pt
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  

  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```

## Interaction: Proportion NU in District (2014)

<br>

```{r, echo=TRUE}
dat$mod = dat$prop_kec_nu
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```



## Interaction: Proportion Muhammadiyah in District (2014)

<br>

```{r, echo=TRUE}
dat$mod = dat$prop_kec_muha
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")

coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)


for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_int_cov(panel_only = TRUE, d = dat)
  
  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s):mod")
  coefs = rbind(coefs,
                data.frame(dv = vec_dvlab[i],
                           coef = out$df$Coefficient[idx],
                           se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

```

<!---BLOCK_LANDSCAPE_STOP--->

# Regression Models with Composite Indices

<br>


```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))


dat$treat = ifelse(dat$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))


```


## Psychometric Properties

<br>

Outgroup rejection:

```{r, echo=FALSE, warning=FALSE, message=FALSE}

### outgroup rejection
cols = c("live_village2pt", "live_neighborhood2pt", "live_room2pt", "worship_house2pt", "intermarriage2pt")

index = which(complete.cases(dat[, cols]))
psych::alpha(dat[index, cols])$total[1]
summary(prcomp(dat[index, cols]))
dat$pca_reject[index] = prcomp(dat[index, cols])$x[,1]
```

<br>

Candidate preferences:

```{r, echo=FALSE, warning=FALSE, message=FALSE}
### candidate preference
cols = c("similar_religion", "similar_ethnicity")

index = which(complete.cases(dat[, cols]))
psych::alpha(dat[index, cols])$total[1]
summary(prcomp(dat[index, cols]))
dat$pca_candidate[index] = prcomp(dat[index, cols])$x[,1]
```

<br>

Ingroup trust and good neighborliness:

```{r, echo=FALSE, warning=FALSE, message=FALSE}
### trust and neighborliness
cols = c("trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

index = which(complete.cases(dat[, cols]))
psych::alpha(dat[index, cols])$total[1]
summary(prcomp(dat[index, cols]))
dat$pca_group[index] = prcomp(dat[index, cols])$x[,1]

```


## Regression

<br>


```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("pca_reject", "pca_candidate", "pca_group")

vec_dvlab = c("Outgroup Rejection",
              "Candidate Preference", 
              "Group Attachment")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  
  
  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_robust_placebo/coefs_composite.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")

# coefs
```

\newpage

<!---BLOCK_LANDSCAPE_START--->

# Treatment Period Variation

<br>

Given that IFLS-4 was fielded in 2007 and IFLS-5 in 2014, the main regression models define treatment as having at least one new mosque constructed between 2008 and 2013, inclusive. The "treatment period variation" analysis adjusts this definition of the treatment period. Four additional models examine treatment effects when the treatment is defined as having at least one new mosque constructed between 2009 and 2013, 2010 and 2013, 2011 and 2013, and 2012 and 2013.

<br>

Tables below show that the results largely hold up even with these different definitions of treatment period.

<br>

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))




```


## Treatment: New Mosques Built from 2009 to 2013, Inclusive

<br>

```{r, echo=FALSE}
dat$treat = ifelse(dat$change_mo_08_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```

Number of districts by treatment status:

```{r, echo=FALSE}

d = dat[!duplicated(dat[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(d$treat))
tab2 = data.frame(round(prop.table(table(d$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status:

```{r, echo=FALSE}

tab1 = data.frame(table(dat$treat) / 2)
tab2 = data.frame(round(prop.table(table(dat$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")
```

<br>

```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_parallel/coefs_09.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")


```

## Treatment: New Mosques Built from 2010 to 2013, Inclusive

<br>

```{r, echo=FALSE}
dat$treat = ifelse(dat$change_mo_09_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```

<br>

Number of districts by treatment status:

```{r, echo=FALSE}

d = dat[!duplicated(dat[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(d$treat))
tab2 = data.frame(round(prop.table(table(d$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status:

```{r, echo=FALSE}

tab1 = data.frame(table(dat$treat) / 2)
tab2 = data.frame(round(prop.table(table(dat$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")
```

<br>


```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_parallel/coefs_10.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")


```


## Treatment: New Mosques Built from 2011 to 2013, Inclusive

<br>

```{r, echo=FALSE}
dat$treat = ifelse(dat$change_mo_10_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```

<br>

Number of districts by treatment status:

```{r, echo=FALSE}

d = dat[!duplicated(dat[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(d$treat))
tab2 = data.frame(round(prop.table(table(d$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status:

```{r, echo=FALSE}

tab1 = data.frame(table(dat$treat) / 2)
tab2 = data.frame(round(prop.table(table(dat$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")
```

<br>


```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_parallel/coefs_11.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")


```


## Treatment: New Mosques Built from 2012 to 2013, Inclusive

<br>

```{r, echo=FALSE}
dat$treat = ifelse(dat$change_mo_11_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

```

<br>

Number of districts by treatment status:

```{r, echo=FALSE}

d = dat[!duplicated(dat[, c("kecid", "treat")]), c("kecid", "treat")]

tab1 = data.frame(table(d$treat))
tab2 = data.frame(round(prop.table(table(d$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")

```

<br>

Number of respondents by treatment status:

```{r, echo=FALSE}

tab1 = data.frame(table(dat$treat) / 2)
tab2 = data.frame(round(prop.table(table(dat$treat)),3))

tab = cbind(tab1, tab2[,2])
colnames(tab) = c("Treatment", "Frequency", "Proportion")

out <- flextable(tab)
out = fontsize(out, size = 7, part = "all")
set_table_properties(out, width = .6, layout = "autofit")
```

<br>


```{r, echo=FALSE, warning=FALSE, message=FALSE}

vec_dv = c("live_village2pt", "live_neighborhood2pt", "live_room2pt",
           "worship_house2pt", "intermarriage2pt", 
           "similar_religion", "similar_ethnicity", 
           "trust_religion2pt", "trust_ethnicity2pt", "helping2pt")

vec_dvlab = c("Village",
              "Neighborhood", 
              "House",
              "House_of_Worship",
              "Intermarriage",
              "Prefer_Co_Religionist", 
              "Prefer_Co_Ethnic",
              "Trust_Muslims", 
              "Trust_Co_Ethnic", 
              "Helping")
                    
coefs = data.frame(dv = NA,
                   coef = NA,
                   se = NA,
                   p = NA)

  
for (i in 1:length(vec_dv)) {
  
  dat$dv = dat[, vec_dv[i]]
  out = regs_basic(panel_only = TRUE, d = dat)

  idx = which(out$df$Parameter=="year2014:treatNew Mosque(s)")
  coefs = rbind(coefs,
              data.frame(dv = vec_dvlab[i],
                   coef = out$df$Coefficient[idx],
                   se = out$df$SE[idx], p = out$df$p[idx]))
  
  out$df$Parameter = change_coefs(out$df$Parameter)
  
  out$df$Coefficient = round(out$df$Coefficient,3)
  out$df$SE = round(out$df$SE,2)
  
  t = data.frame(out$df)
  t$SE = paste0("\n(",t$SE,")")
  
  t$Coefficient[which(t$p <= .001)] = paste0(t$Coefficient[which(t$p <= .001)], "*")
  t$Coefficient[which(t$p <= .01)] = paste0(t$Coefficient[which(t$p <= .01)], "*")
  t$Coefficient[which(t$p <= .05)] = paste0(t$Coefficient[which(t$p <= .05)], "*")
  
  t$Coefficient = paste0(t$Coefficient, t$SE)  

  t[nrow(t),2] = gsub("\\n\\(NA\\)", "", t[nrow(t),2])
  t[nrow(t)-1,2] = gsub("\\n\\(NA\\)", "", t[nrow(t)-1,2])
  
  if (i == 1) {
    tab = t[,c(1,2)]
  } else {
    tab = cbind(tab,t[,2])
  }
  
}

colnames(tab)[2:ncol(tab)] = vec_dvlab

coefs = coefs[-1,]
write.csv(coefs, file = "./coefs_parallel/coefs_12.csv", row.names = FALSE)


tab = flextable(data.frame(tab))
tab = add_footer_lines(tab, "*** p<.001 ** p<.01 * p<.05")
tab = fontsize(tab, size = 7, part = "all")
tab = align(tab, j = 2:(length(vec_dv)+1), align = "center", part = "body")
set_table_properties(tab, width = 1, layout = "autofit")


```


<!---BLOCK_LANDSCAPE_STOP--->

# Parallel Trends (All Violent Incidents)

<br>

One of the assumptions of the Difference-in-Differences (DiD) approach is that the control and treated units exhibit parallel trends prior to the application of the treatment (in this case, prior to 2008). Since IFLS waves prior to Wave 4 did not include questions related to social and political attitudes, it is therefore impossible to explicitly examine this assumption using IFLS data.

<br>

However, the National Violence Monitoring System (NVMS), developed jointly by the World Bank and the Indonesian Government, provides an alternative means to check this assumption to some extent. Using this dataset, we can examine whether the treated districts were disproportionately susceptible to social conflicts and violence compared to the control districts.

<br>

The figures below present the average number of violent incidents in the control and treated districts between 1998 and 2014, using different operational definitions of the treatment period. As the previous section demonstrated that these different operationalizations yield practically identical substantive results, the primary interest here is to determine whether the trends for the treated and control districts are reasonably parallel prior to the onset of the treatment. The figures show that this expectation is met in all cases, as the lines are practically overlapping.

<br>


```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

```

```{r, echo=FALSE, message=FALSE, warning=FALSE}

## collapse by kecamatan
cols = c("kecid14", "kecid", 
         "change_mo_07_13", "change_mo_08_13", "change_mo_09_13", "change_mo_10_13", 
         "change_mo_11_13", "change_mo_12_13")
kec = dat[, cols]
kec = unique.data.frame(kec)

## prepare nvms
load("data/nvms_kec.RData")


## crosswalk and merge
cw_nvms = read.xlsx2("data/crosswalk_kec_bps2014_nvms.xlsx", 
                     sheetIndex=1)

cw_nvms = subset(cw_nvms, subset = src_nvms != "")

temp = merge(x = kec, y = cw_nvms,
             all.x = TRUE, all.y = FALSE)

temp = temp[!duplicated(temp$kecid),]

kec = merge(x = temp, y = nvms_kec,
            all = FALSE)

```

\newpage

## Treatment: New Mosques (2008-2013, Inclusive)

<br>

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6.5}

kec$Treatment = ifelse(kec$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
kec$Treatment = factor(kec$Treatment, levels = c("No New Mosque", "New Mosque(s)"))
# table(kec$Treatment)

file = "./coefs_robust_placebo/coefs_main_models.csv"

p <- plot_parallel(start = 2008, file)

p$p_all
p$p_eff

```

\newpage

## Treatment: New Mosques (2009-2013, Inclusive)

<br>

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6.5}

kec$Treatment = ifelse(kec$change_mo_08_13 > 0, "New Mosque(s)", "No New Mosque")
kec$Treatment = factor(kec$Treatment, levels = c("No New Mosque", "New Mosque(s)"))
# table(kec$Treatment)

file = "./coefs_parallel/coefs_09.csv"

p <- plot_parallel(start = 2009, file)

p$p_all
p$p_eff

```

\newpage

## Treatment: New Mosques (2010-2013, Inclusive)

<br>

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6.5}

kec$Treatment = ifelse(kec$change_mo_09_13 > 0, "New Mosque(s)", "No New Mosque")
kec$Treatment = factor(kec$Treatment, levels = c("No New Mosque", "New Mosque(s)"))
# table(kec$Treatment)

file = "./coefs_parallel/coefs_10.csv"

p <- plot_parallel(start = 2010, file)

p$p_all
p$p_eff

```

\newpage

## Treatment: New Mosques (2011-2013, Inclusive)

<br>

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6.5}

kec$Treatment = ifelse(kec$change_mo_10_13 > 0, "New Mosque(s)", "No New Mosque")
kec$Treatment = factor(kec$Treatment, levels = c("No New Mosque", "New Mosque(s)"))
# table(kec$Treatment)

file = "./coefs_parallel/coefs_11.csv"

p <- plot_parallel(start = 2011, file)

p$p_all
p$p_eff

```

\newpage

## Treatment: New Mosques (2012-2013, Inclusive)

<br>

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6.5}

kec$Treatment = ifelse(kec$change_mo_11_13 > 0, "New Mosque(s)", "No New Mosque")
kec$Treatment = factor(kec$Treatment, levels = c("No New Mosque", "New Mosque(s)"))
# table(kec$Treatment)

file = "./coefs_parallel/coefs_12.csv"

p <- plot_parallel(start = 2012, file)

p$p_all
p$p_eff

```

\newpage

# Parallel Trends (Religion-Related Violent Incidents)

<br>

The figures below present the average numbers of religion-related violent incidents in the control and treated districts between 1998 and 2014, using different operational definitions of treatment period. The trends in the control and treated districts appear reasonably parallel when the treatment is defined as at least one new mosque constructed between 2008–2013 and 2009–2013. The trends are less parallel for the other treatment periods. However, the patterns contradict the documented treatment effects, as during these periods, the control districts experienced more religion-related violent incidents.



```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.rData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

```

```{r, echo=FALSE, message=FALSE, warning=FALSE}

## collapse by kecamatan
cols = c("kecid14", "kecid", 
         "change_mo_07_13", "change_mo_08_13", "change_mo_09_13", "change_mo_10_13", 
         "change_mo_11_13", "change_mo_12_13")
kec = dat[, cols]
kec = unique.data.frame(kec)

## prepare nvms
load("data/nvms_kec.RData")


## crosswalk and merge
cw_nvms = read.xlsx2("data/crosswalk_kec_bps2014_nvms.xlsx", 
                     sheetIndex=1)

cw_nvms = subset(cw_nvms, subset = src_nvms != "")

temp = merge(x = kec, y = cw_nvms,
             all.x = TRUE, all.y = FALSE)

temp = temp[!duplicated(temp$kecid),]

kec = merge(x = temp, y = nvms_kec,
            all = FALSE)

```

\newpage

## Treatment: New Mosques (2008-2013, Inclusive)

<br>

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6.5}

kec$Treatment = ifelse(kec$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
kec$Treatment = factor(kec$Treatment, levels = c("No New Mosque", "New Mosque(s)"))
# table(kec$Treatment)

file = "./coefs_robust_placebo/coefs_main_models.csv"

p <- plot_parallel(start = 2008, file)

p$p_relmean
# p$p2

```
\newpage

## Treatment: New Mosques (2009-2013, Inclusive)

<br>

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6.5}

kec$Treatment = ifelse(kec$change_mo_08_13 > 0, "New Mosque(s)", "No New Mosque")
kec$Treatment = factor(kec$Treatment, levels = c("No New Mosque", "New Mosque(s)"))
# table(kec$Treatment)

file = "./coefs_parallel/coefs_09.csv"

p <- plot_parallel(start = 2009, file)

p$p_relmean
# p$p2

```

\newpage

## Treatment: New Mosques (2010-2013, Inclusive)

<br>

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6.5}

kec$Treatment = ifelse(kec$change_mo_09_13 > 0, "New Mosque(s)", "No New Mosque")
kec$Treatment = factor(kec$Treatment, levels = c("No New Mosque", "New Mosque(s)"))
# table(kec$Treatment)

file = "./coefs_parallel/coefs_10.csv"

p <- plot_parallel(start = 2010, file)

p$p_relmean
# p$p2

```

\newpage

## Treatment: New Mosques (2011-2013, Inclusive)

<br>

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6.5}

kec$Treatment = ifelse(kec$change_mo_10_13 > 0, "New Mosque(s)", "No New Mosque")
kec$Treatment = factor(kec$Treatment, levels = c("No New Mosque", "New Mosque(s)"))
# table(kec$Treatment)

file = "./coefs_parallel/coefs_11.csv"

p <- plot_parallel(start = 2011, file)

p$p_relmean
# p$p2

```

\newpage

## Treatment: New Mosques (2012-2013, Inclusive)

<br>

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=6.5}

kec$Treatment = ifelse(kec$change_mo_11_13 > 0, "New Mosque(s)", "No New Mosque")
kec$Treatment = factor(kec$Treatment, levels = c("No New Mosque", "New Mosque(s)"))
# table(kec$Treatment)

file = "./coefs_parallel/coefs_12.csv"

p <- plot_parallel(start = 2012, file)

p$p_relmean
# p$p2

```

